% optimize the network based on sparseness measure
function optimize_map

global Ne IN N A AI Achand best

% call simulator once to define the area sizes as globals.
spnet_som( );


i=1;


radius=length(A)/2
conn_data(i,:)=[0.5*radius 0 5 radius 3 5 .05 0 10];i=i+1;

radius=length(AI)/2;
conn_data(i,:)=[0.5*radius 0 5 radius 3 5 .05 0 15];i=i+1;

radius=length(AI)/4;
conn_data(i,:)=[0.5*radius 0 2 radius 1 2 .05 0 15];i=i+1;

radius=length(Achand)/4;
conn_data(i,:)=[radius 0 2 radius 1 2 .05 0 10];i=i+1;

radius = length(A)/6;
conn_data(i,:)=[0.5*radius 0  1 radius 1 2 0.0 0 6];i=i+1;

radius = length(A)/6;
conn_data(i,:)=[0.5*radius 0 -6 radius 1 1 0.0 -6 0];i=i+1;
%connect(AI,A,'gauss', 'std', conn(1), 'mininitw', conn(2), 'maxinitw', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delaymax', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxw', conn(9));

lowest_dist = 10^11;
lowest(1) = -1;
lowest(2) = -1;

bad_count = 0;


target_mfr = zeros(size(A));
target_mfr(24-4:40+4)=(gaussian_1d(4,-12:12))*300;


for i = 1:100000
    
    stop_sec=1;
    start_recording_sec= max(0,stop_sec - 20);
    sdat = spnet_som(conn_data);
    save_seconds = stop_sec - start_recording_sec ;
    
    % bin the spike data by pattern.
    
    S=sparse( sdat(:,2), sdat(:,1), 1, N , stop_sec*1000 );

	% Gather mean rates for stats.
    epoch_ms = 1000;
    patterns = 1;
    isi_ms = epoch_ms/patterns;
    
    t=start_recording_sec*1000 + 1;
    count = zeros(1,patterns);
    mfr = zeros(patterns,N);
    for i = start_recording_sec*1000:epoch_ms:stop_sec*1000-1
        for p = 1:patterns
            mfr(p, : ) = mfr(p, :) + (1000/isi_ms)*full(sum( S( :, t+20:(t+isi_ms-1)), 2))';
            count(p)=count(p)+1;
            t=t+isi_ms;
        end
    end
    mfr=mfr./(repmat( count', 1, size(mfr,2) ) + .000001);
    
    
    % Look at sparsness measures
    patterns = 1;

%     
%     j= 1;
%     for i = max(IN)+1:Ne
%         lifetime_sparseness(j) = sparseness( mfr(:,i), 40 ); % 40 Hz is max
%         j= j + 1;
%     end
%     disp(['Mean lifetime sparseness: ' num2str( mean(lifetime_sparseness))]);
%     
%     for i = 1:patterns
%         population_sparseness(i) = sparseness(mfr(i,max(IN)+1:Ne), 40);
%     end
%     disp(['Mean population sparseness: ' num2str( mean(population_sparseness))]);
        
    
    
    % save the best
%    d= sqrt( (mean(population_sparseness)-0.5)^2 )
%d= sqrt( (mean(population_sparseness)-0.5)^2 + (mean(lifetime_sparseness)-0.5)^2)

    d = sqrt( sum((mfr(:,max(IN)+1:Ne) - target_mfr ).^2 ))

    if d < lowest_dist
        bad_count = 0;
        lowest_dist = d;
        lowest = mfr(:,max(IN)+1:Ne);
        %lowest(2) = mean(lifetime_sparseness);
        best_params = conn_data;
        best_mfr = mfr;
        save 'best' lowest best_params best_mfr;
    else
        bad_count = bad_count + 1;
        
        if bad_count > 5
            conn_data = best_params; % start again from best spot.
        end
    end
        
    
    
    % Tweek the parameters 
    
    conn_data = tweak( conn_data );

    lowest
    conn_data
end



function new_conn = tweak(conn_data)

global Ne IN N A AI



new_conn(1,:) = conn_data(1,:) + [ rand*5-2.5 0 rand*2-1 rand*2-1 round(rand*2-1) round(rand*2-1) rand*.01-.005 rand*.1 rand*1];
new_conn(2,:) = conn_data(2,:) + [ rand*5-2.5 0 rand*2-1 rand*2-1 round(rand*2-1) round(rand*2-1) rand*.01-.005 rand*.1 rand*1];
new_conn(3,:) = conn_data(3,:) + [ rand*5-2.5 0 rand*2-1 rand*2-1 round(rand*2-1) round(rand*2-1) rand*.01-.005 rand*.1 rand*1];
new_conn(4,:) = conn_data(4,:) + [ rand*5-2.5 0 rand*2-1 rand*2-1 round(rand*2-1) round(rand*2-1) rand*.01-.005 rand*.1 rand*1];
new_conn(5,:) = conn_data(5,:) + [ rand*5-2.5 0 rand*2-1 rand*2-1 round(rand*2-1) round(rand*2-1) rand*.01-.005 rand*.1 rand*1];
new_conn(6,:) = conn_data(6,:) + [ rand*5-2.5 0 rand*2-1 rand*2-1 round(rand*2-1) round(rand*2-1) rand*.01-.005 rand*.1 rand*1];

% max wt should be >= min wt.
bad_indices = conn_data(:,9) < new_conn(:,8);
new_conn( bad_indices, 8 ) = new_conn(bad_indices, 9); 

bad_indices = conn_data(:,6) < new_conn(:,5);
new_conn( bad_indices, 5 ) = new_conn(bad_indices, 6); 


max_conn_data(1,:) = [30 5 20 length(A)/2 5 5 .1 5 20];
max_conn_data(2,:) = [30 5 20 length(AI)/2 5 5 .1 5 20];
max_conn_data(3,:) = [30 5 20 length(AI)/2 5 5 .1 5 20];
max_conn_data(4,:) = [30 5 20 length(A)/2 5 5 .1 5 20];
max_conn_data(5,:) = [30 5 20 length(A)/2 5 5 .1 5 20];
max_conn_data(6,:) = [30 0 0 length(A)/2 5 5 .1 5  0];

min_conn_data(1,:) = [.5 0 1 1 1 1 0 0 0];
min_conn_data(2,:) = [.5 0 1 1 1 1 0 0 0];
min_conn_data(3,:) = [.5 0 1 1 1 1 0 0 0];
min_conn_data(4,:) = [.5 0 1 1 1 1 0 0 0];
min_conn_data(5,:) = [.5 0 1 1 1 1 0 0 0];
min_conn_data(6,:) = [.5 -15 -20 1 1 1 0 -20 0];

new_conn = max(min_conn_data,min(max_conn_data,new_conn));
%connect(IN,A,'gauss', 'std', conn(1), 'mininitw', conn(2), 'maxinitw', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delaymax', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxw', conn(9));



